{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 5: Introduction to Probability Theory (Part III)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Objectives\n",
    "\n",
    "+ To introduce continuous random variables.\n",
    "+ To introduce the cumulative distribution function and its properties.\n",
    "+ To introduce the probability density function and its properties.\n",
    "+ To introduce the expectation of continuous random variables.\n",
    "+ To introduce the concept of the joint probability density function.\n",
    "+ To present some simple analytical examples Bayesian inference."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Readings\n",
    "\n",
    "+ These notes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set_context('talk')\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Something more on probability spaces\n",
    "\n",
    "It turns out that when $\\Omega$ is a continuous space, like $\\mathbb{R}$ for example, it is not possible to take $\\mathcal{F}$ to be all the subsets of $\\Omega$ (you need to take a measure-theoretic probability theory class to understand why this is).\n",
    "However, many \"nice\" subsets of $\\Omega$ are usually in $\\mathcal{F}$.\n",
    "For example, in the case of $\\Omega=\\mathbb{R}$, $\\mathcal{F}$ can include all intervals, and any countable unions and intersections of intervals.\n",
    "That's a lot of sets.\n",
    "\n",
    "In any case, $\\mathcal{F}$ must satisfy certain properties for everything to be well-defined.\n",
    "These properties are:\n",
    "+ $\\Omega \\in \\mathcal{F}$\n",
    "+ For any $A$ in $\\mathcal{F}$, the complement $A^c$ is in $\\mathcal{F}$.\n",
    "+ For any $A_1,A_2,\\dots$ in $\\mathcal{F}$, the union $\\cup_n A_n$ is in $\\mathcal{F}$.\n",
    "When a set of subsets $\\mathcal{F}$ satisfies these properties, we say that it forms a $\\sigma$-algebra."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Continuous random variables\n",
    "\n",
    "Consider a probability space $(\\Omega, \\mathcal{F}, \\mathbb{P})$.\n",
    "A continuous random variable models the result of an experiment that can potentially take infinitely many values.\n",
    "That is, it is a function\n",
    "$$\n",
    "X:\\Omega \\rightarrow \\mathbb{R}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The cumulative distribution function\n",
    "\n",
    "Let $X$ be a continuous random variable.\n",
    "Its cumulative distribution function (CDF) $F_X(x)$ gives the probability that $X$ is smaller than $x$:\n",
    "$$\n",
    "F_X(x) := \\mathbb{P}(X\\le x) = \\mathbb{P}\\left(\\left\\{\\omega: X(\\omega) \\le x\\right\\}\\right).\n",
    "$$\n",
    "\n",
    "### Properties of the cumulative distribution function\n",
    "\n",
    "+ $F_X(x)$ is an increasing function.\n",
    "+ $F_X(-\\infty) = 0$.\n",
    "+ $F_X(+\\infty) = 1$.\n",
    "+ $\\mathbb{P}(a\\le X \\le b) = F_X(b) - F_X(a)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The probability density function\n",
    "\n",
    "The probability density function (PDF) is a \"function\" $f_X(x)$ that can give us the probability that $X$ is in any \"good\" subset $A$ of $\\mathbb{R}$ as follows:\n",
    "$$\n",
    "\\mathbb{P}(X\\in A) = \\int_A f_X(x) dx.\n",
    "$$\n",
    "Note that certain random variables may not have a PDF that is a function.\n",
    "That's why I put the word \"function\" in quotes.\n",
    "However, if you allow the PDF to include Dirac's $\\delta$, then any random variable has a PDF.\n",
    "We will ignore this complication for the moment.\n",
    "\n",
    "In this class, we will simplify the notation and we will be writing:\n",
    "$$\n",
    "p(x) \\equiv f_X(x),\n",
    "$$\n",
    "when there is no ambiguity.\n",
    "\n",
    "### Properties of the probability density function\n",
    "+ $p(x) \\ge 0$ for all $x$.\n",
    "+ $\\int_{-\\infty}^{\\infty} p(x) dx = 1$.\n",
    "+ The derivative of the CDF is the PDF, i.e., $F_X'(x) = p(x)$.\n",
    "\n",
    "### Dirac's delta and a unified view of all random variables\n",
    "\n",
    "Dirac's $\\delta$ is a special function, actually called a distribution, which is defined as follows:\n",
    "$$\n",
    "\\delta(x) = 0,\n",
    "$$\n",
    "for $x\\not=0$, and\n",
    "$$\n",
    "\\int_{-\\infty}^\\infty \\delta(x) dx = 0\n",
    "$$\n",
    "You can think of Dirac's $\\delta$ as the PDF of a discrete random variable that takes the value $0$ with probability one.\n",
    "\n",
    "Using Dirac's $\\delta$ you can make any discrete random variable you want look like a continuous random variable.\n",
    "For example, consider a Categorical random variable taking values $x_1,\\dots,x_K$ with probabilities $p_1,\\dots,p_K$.\n",
    "The PDF of this random variable can be written as:\n",
    "$$\n",
    "p(x) = \\sum_{k=1}^Kp_k\\delta(x - x_k).\n",
    "$$\n",
    "\n",
    "The most general random variable one can think of has a PDF that consists of two parts:\n",
    "$$\n",
    "p(x) = f^n_X(x) + f^{\\delta}_X(x),\n",
    "$$\n",
    "a part $f^n_X(x)$ that is a nice proper function and a part $f^{\\delta}_X(x)$ that consists of a weighted sum of Dirac $\\delta$'s. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Expectations of continuous random variables\n",
    "\n",
    "The expectation of a continuous random variable is:\n",
    "$$\n",
    "\\mathbb{E}[X] = \\int_{-\\infty}^\\infty x p(x)dx.\n",
    "$$\n",
    "Its properties are the same as the expectation of a discrete random variable.\n",
    "\n",
    "The expectation of a function of the random variable is:\n",
    "$$\n",
    "\\mathbb{E}[f(X)] = \\int_{-\\infty}^\\infty f(x)p(x)dx.\n",
    "$$\n",
    "\n",
    "In general, using the Dirac $\\delta$ trick we introduced above, there is no need to differentiate between integrations and summations.\n",
    "We can use integration for both continuous and discrete random variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: The uniform distribution\n",
    "\n",
    "The uniform distribution is the most common continuous distribution.\n",
    "It corresponds to a random variable that is equally likely to take a value within a given interval.\n",
    "We write:\n",
    "$$\n",
    "X\\sim U([0,1]),\n",
    "$$\n",
    "and we read $X$ follows a uniform distribution taking values in $[0,1]$.\n",
    "\n",
    "The probability density of the uniform is constant in $[0,1]$ and zero outside it.\n",
    "We have:\n",
    "$$\n",
    "p(x) := U(x|[0,1]) := f_X(x) = \\begin{cases}\n",
    "1,&\\;0\\le x \\le 1,\\\\\n",
    "0,&\\;\\text{otherwise}.\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "The cumulative distribution funciton of the uniform is:\n",
    "$$\n",
    "F_X(x) = \\mathbb{P}(X \\le x) = \\int_0^x f_X(u) du = \\int_0^x du = x.\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbb{P}(a \\le X \\le b) = F_X(b) - F_X(a) = b - a.\n",
    "$$\n",
    "\n",
    "The expectation of the uniform is:\n",
    "$$\n",
    "\\mathbb{E}[X] = \\int_0^1 xdx = \\frac{1}{2}.\n",
    "$$\n",
    "\n",
    "The variance of the uniform is:\n",
    "$$\n",
    "\\mathbb{V}[X] = \\mathbb{E}[X^2] - \\left(\\mathbb{E}[X]\\right)^2 = \\frac{1}{3} - \\frac{1}{4} = \\frac{1}{12}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.97837346, 0.625879  , 0.52290077, 0.51240399, 0.04980333,\n",
       "       0.59070721, 0.99822886, 0.88110047, 0.77926415, 0.80169289,\n",
       "       0.63377297, 0.4330944 , 0.39263228, 0.87655016, 0.71793991,\n",
       "       0.96914239, 0.98791565, 0.03918648, 0.24459545, 0.32064884,\n",
       "       0.07865616, 0.6215369 , 0.44383284, 0.98624492, 0.23687563,\n",
       "       0.76983114, 0.76920715, 0.12850447, 0.5412241 , 0.54202371,\n",
       "       0.56063597, 0.46389922, 0.41674185, 0.03591847, 0.09757618,\n",
       "       0.15020877, 0.62910003, 0.95689217, 0.59937696, 0.03754368,\n",
       "       0.45902144, 0.18677698, 0.91152939, 0.84733429, 0.57136117,\n",
       "       0.38750976, 0.75984272, 0.82609542, 0.95836646, 0.86262593,\n",
       "       0.71400927, 0.35642015, 0.20075114, 0.23625662, 0.50861275,\n",
       "       0.92449095, 0.59173485, 0.21624696, 0.23393397, 0.76204738,\n",
       "       0.10804259, 0.27232046, 0.51790518, 0.96678481, 0.43760293,\n",
       "       0.56561303, 0.95582303, 0.75536118, 0.40260221, 0.83693531,\n",
       "       0.78007319, 0.28611426, 0.64475304, 0.70836801, 0.23042198,\n",
       "       0.91461071, 0.54027546, 0.20290475, 0.10600445, 0.36341161,\n",
       "       0.58776198, 0.96599455, 0.4787342 , 0.89632266, 0.11829211,\n",
       "       0.53011026, 0.10231121, 0.48088843, 0.68610058, 0.05500749,\n",
       "       0.0664882 , 0.83686234, 0.87464224, 0.5527214 , 0.98598549,\n",
       "       0.83999942, 0.50744833, 0.78022083, 0.94580147, 0.66221155])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's demonstrate how you can sample from the uniform\n",
    "import scipy.stats as st\n",
    "X = st.uniform()\n",
    "X.rvs(size=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.70121505, 0.91085994, 0.69158752, 0.11946638, 0.09728609,\n",
       "       0.49065546, 0.27969302, 0.21582279, 0.93441454, 0.41349678,\n",
       "       0.47870551, 0.73636654, 0.13051002, 0.37208425, 0.3116429 ,\n",
       "       0.82258973, 0.23919661, 0.22838676, 0.5542224 , 0.53016444,\n",
       "       0.35530159, 0.09605987, 0.5549514 , 0.60081403, 0.73167165,\n",
       "       0.78542101, 0.32674251, 0.03529146, 0.59239157, 0.90631834,\n",
       "       0.61845217, 0.60454155, 0.17810513, 0.62028272, 0.15576036,\n",
       "       0.65303787, 0.7683246 , 0.49567772, 0.54411001, 0.22508122,\n",
       "       0.84317487, 0.8576526 , 0.66822961, 0.41742108, 0.00306515,\n",
       "       0.87628444, 0.02520404, 0.98943986, 0.00223345, 0.17440434,\n",
       "       0.52514865, 0.00128126, 0.92851971, 0.00146251, 0.49170257,\n",
       "       0.75837021, 0.88191865, 0.29204905, 0.20065024, 0.07695563,\n",
       "       0.1160243 , 0.55381338, 0.82336444, 0.0619979 , 0.49091839,\n",
       "       0.27736099, 0.35720052, 0.12455454, 0.39670171, 0.98057418,\n",
       "       0.65956781, 0.22232854, 0.44287494, 0.4036119 , 0.69792346,\n",
       "       0.10159559, 0.69955144, 0.50623539, 0.7209267 , 0.69582729,\n",
       "       0.16009671, 0.69380665, 0.61887288, 0.66828418, 0.54099226,\n",
       "       0.52855516, 0.62526337, 0.08661942, 0.02498284, 0.9058553 ,\n",
       "       0.43726005, 0.85914546, 0.3964603 , 0.05790563, 0.93392332,\n",
       "       0.51886644, 0.2245634 , 0.03197055, 0.00429038, 0.06374447])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# An alternative way is to use the functionality of numpy\n",
    "np.random.rand(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions\n",
    "\n",
    "+ Modify the code above so that you sample from a uniform distribution taking values in $[2, 5]$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Joint probability density function\n",
    "\n",
    "Let $X$ and $Y$ be two random variables.\n",
    "There joint probability density $f_{X,Y}(x,y)$ is the function that can give us the probability that the pair $(X,Y)$ belongs to any \"good\" subset $A$ of $\\mathbb{R}^2$ as follows:\n",
    "$$\n",
    "\\mathbb{P}\\left((X,Y)\\in A\\right) = \\int\\int_{A} f_{X,Y}(x,y)dxdy.\n",
    "$$\n",
    "Of course, we will be writing:\n",
    "$$\n",
    "p(x,y) := f_{X,Y}(x,y),\n",
    "$$\n",
    "when there is no ambiguity.\n",
    "\n",
    "If you integrate one of the variables out of the joint, you get the PDF of the other variable.\n",
    "For example:\n",
    "$$\n",
    "p(x) = \\int_{-\\infty}^\\infty p(x,y) dy,\n",
    "$$\n",
    "and\n",
    "$$\n",
    "p(y) = \\int_{-\\infty}^\\infty p(x, y) dx.\n",
    "$$\n",
    "\n",
    "For many random variables $X_1,\\dots,X_N$ their joint PDF $p(x_1,\\dot,x_N)$ is similarly defined.\n",
    "Again, integrating out some variables gives the pdf of the others."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conditioning a random variable on another\n",
    "\n",
    "Consider two random variables $X$ and $Y$.\n",
    "If we had observed that $Y=y$, how would this change the PDF of $X$?\n",
    "The answer is given via Bayes' rule.\n",
    "The PDF of $X$ conditioned on $Y=y$ is:\n",
    "$$\n",
    "p(x|y) = \\frac{p(x,y)}{p(y)}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Inferring the probability of a coin toss from data (1/3)\n",
    "\n",
    "This is our first Bayesian inference example!\n",
    "\n",
    "Let's say that we run a coin toss experiment $N$ times and we wish to figure out the probability of heads.\n",
    "We will bring into the picture all the mathematical machinery we have introduced so far.\n",
    "Let's say that the data we have observe are $x_1,\\dots,x_N$.\n",
    "For notational convenience we will be writing:\n",
    "$$\n",
    "x_{1:N} := (x_1,\\dots,x_N).\n",
    "$$\n",
    "\n",
    "First, let's start with the probability of success of the coin toss.\n",
    "Let's call it $\\theta$.\n",
    "How can we describe our uncertainty about it?\n",
    "We have to assign a *prior* probability distribution on it.\n",
    "Let's say that we don't know anything about it except that it must be between 0 and 1.\n",
    "What distribution should we assign?\n",
    "Of course, a uniform distribution:\n",
    "$$\n",
    "\\theta \\sim U([0,1]).\n",
    "$$\n",
    "Second, each coin toss experiment corresponds to an independent Bernoulli variable with the same probability of success $\\theta$.\n",
    "We write:\n",
    "$$\n",
    "X_n | \\theta \\sim \\operatorname{Bernoulli}(\\theta),\n",
    "$$\n",
    "for $n=1,\\dots,N$.\n",
    "Note that these random variables depend on $\\theta$.\n",
    "That's why we are conditioning like this.\n",
    "\n",
    "Before proceeding with the mathematics, let's draw the graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: coin_toss_bayes_1 Pages: 1 -->\n",
       "<svg width=\"62pt\" height=\"116pt\"\n",
       " viewBox=\"0.00 0.00 62.00 116.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 112)\">\n",
       "<title>coin_toss_bayes_1</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-112 58,-112 58,4 -4,4\"/>\n",
       "<!-- theta -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>theta</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"27\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"23\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">θ</text>\n",
       "</g>\n",
       "<!-- X1 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>X1</title>\n",
       "<ellipse fill=\"#d3d3d3\" stroke=\"#000000\" cx=\"27\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"18.5\" y=\"-15.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">X</text>\n",
       "<text text-anchor=\"start\" x=\"29.5\" y=\"-15.3\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\" fill=\"#000000\">1</text>\n",
       "</g>\n",
       "<!-- theta&#45;&gt;X1 -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>theta&#45;&gt;X1</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M27,-71.8314C27,-64.131 27,-54.9743 27,-46.4166\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"30.5001,-46.4132 27,-36.4133 23.5001,-46.4133 30.5001,-46.4132\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x105c31410>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# First for one coin toss\n",
    "from graphviz import Digraph\n",
    "gc = Digraph('coin_toss_bayes_1')\n",
    "gc.node('theta', label='<&theta;>')\n",
    "gc.node('X1', label='<X<sub>1</sub>>', style='filled')\n",
    "gc.edge('theta', 'X1')\n",
    "gc.render('coin_toss_bayes_1', format='png')\n",
    "gc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: coin_toss_bayes_2 Pages: 1 -->\n",
       "<svg width=\"134pt\" height=\"116pt\"\n",
       " viewBox=\"0.00 0.00 134.00 116.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 112)\">\n",
       "<title>coin_toss_bayes_2</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-112 130,-112 130,4 -4,4\"/>\n",
       "<!-- theta -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>theta</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"63\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"59\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">θ</text>\n",
       "</g>\n",
       "<!-- X1 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>X1</title>\n",
       "<ellipse fill=\"#d3d3d3\" stroke=\"#000000\" cx=\"27\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"18.5\" y=\"-15.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">X</text>\n",
       "<text text-anchor=\"start\" x=\"29.5\" y=\"-15.3\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\" fill=\"#000000\">1</text>\n",
       "</g>\n",
       "<!-- theta&#45;&gt;X1 -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>theta&#45;&gt;X1</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M54.2854,-72.5708C50.0403,-64.0807 44.8464,-53.6929 40.1337,-44.2674\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"43.237,-42.6477 35.6343,-35.2687 36.976,-45.7782 43.237,-42.6477\"/>\n",
       "</g>\n",
       "<!-- X2 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>X2</title>\n",
       "<ellipse fill=\"#d3d3d3\" stroke=\"#000000\" cx=\"99\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"90.5\" y=\"-15.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">X</text>\n",
       "<text text-anchor=\"start\" x=\"101.5\" y=\"-15.3\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\" fill=\"#000000\">2</text>\n",
       "</g>\n",
       "<!-- theta&#45;&gt;X2 -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>theta&#45;&gt;X2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M71.7146,-72.5708C75.9597,-64.0807 81.1536,-53.6929 85.8663,-44.2674\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"89.024,-45.7782 90.3657,-35.2687 82.763,-42.6477 89.024,-45.7782\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1a18b3ec90>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Now for two coin tosses\n",
    "gc2 = Digraph('coin_toss_bayes_2')\n",
    "gc2.node('theta', label='<&theta;>')\n",
    "gc2.node('X1', label='<X<sub>1</sub>>', style='filled')\n",
    "gc2.node('X2', label='<X<sub>2</sub>>', style='filled')\n",
    "gc2.edge('theta', 'X1')\n",
    "gc2.edge('theta', 'X2')\n",
    "gc2.render('coin_toss_bayes_2', format='png')\n",
    "gc2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: coin_toss_bayes_3 Pages: 1 -->\n",
       "<svg width=\"206pt\" height=\"116pt\"\n",
       " viewBox=\"0.00 0.00 206.00 116.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 112)\">\n",
       "<title>coin_toss_bayes_3</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-112 202,-112 202,4 -4,4\"/>\n",
       "<!-- theta -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>theta</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"99\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"95\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">θ</text>\n",
       "</g>\n",
       "<!-- X1 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>X1</title>\n",
       "<ellipse fill=\"#d3d3d3\" stroke=\"#000000\" cx=\"27\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"18.5\" y=\"-15.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">X</text>\n",
       "<text text-anchor=\"start\" x=\"29.5\" y=\"-15.3\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\" fill=\"#000000\">1</text>\n",
       "</g>\n",
       "<!-- theta&#45;&gt;X1 -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>theta&#45;&gt;X1</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M83.7307,-74.7307C73.803,-64.803 60.6847,-51.6847 49.5637,-40.5637\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"51.7933,-37.8436 42.2473,-33.2473 46.8436,-42.7933 51.7933,-37.8436\"/>\n",
       "</g>\n",
       "<!-- X2 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>X2</title>\n",
       "<ellipse fill=\"#d3d3d3\" stroke=\"#000000\" cx=\"99\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"90.5\" y=\"-15.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">X</text>\n",
       "<text text-anchor=\"start\" x=\"101.5\" y=\"-15.3\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\" fill=\"#000000\">2</text>\n",
       "</g>\n",
       "<!-- theta&#45;&gt;X2 -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>theta&#45;&gt;X2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M99,-71.8314C99,-64.131 99,-54.9743 99,-46.4166\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"102.5001,-46.4132 99,-36.4133 95.5001,-46.4133 102.5001,-46.4132\"/>\n",
       "</g>\n",
       "<!-- X3 -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>X3</title>\n",
       "<ellipse fill=\"#d3d3d3\" stroke=\"#000000\" cx=\"171\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"162.5\" y=\"-15.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">X</text>\n",
       "<text text-anchor=\"start\" x=\"173.5\" y=\"-15.3\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\" fill=\"#000000\">3</text>\n",
       "</g>\n",
       "<!-- theta&#45;&gt;X3 -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>theta&#45;&gt;X3</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M114.2693,-74.7307C124.197,-64.803 137.3153,-51.6847 148.4363,-40.5637\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"151.1564,-42.7933 155.7527,-33.2473 146.2067,-37.8436 151.1564,-42.7933\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1a18b4db90>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Now for three coin tosses\n",
    "gc3 = Digraph('coin_toss_bayes_3')\n",
    "gc3.node('theta', label='<&theta;>')\n",
    "gc3.node('X1', label='<X<sub>1</sub>>', style='filled')\n",
    "gc3.node('X2', label='<X<sub>2</sub>>', style='filled')\n",
    "gc3.node('X3', label='<X<sub>3</sub>>', style='filled')\n",
    "gc3.edge('theta', 'X1')\n",
    "gc3.edge('theta', 'X2')\n",
    "gc3.edge('theta', 'X3')\n",
    "gc3.render('coin_toss_bayes_3', format='png')\n",
    "gc3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because this notation gets a little bit too involved, we introduce the so called [plate notation](https://en.wikipedia.org/wiki/Plate_notation).\n",
    "Whategver is inside the subgrpah indicated by the box is supposed to be repeated as many times as indicated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: coin_toss_bayes_plate Pages: 1 -->\n",
       "<svg width=\"94pt\" height=\"155pt\"\n",
       " viewBox=\"0.00 0.00 94.00 155.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 151)\">\n",
       "<title>coin_toss_bayes_plate</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-151 90,-151 90,4 -4,4\"/>\n",
       "<g id=\"clust1\" class=\"cluster\">\n",
       "<title>cluster_0</title>\n",
       "<polygon fill=\"none\" stroke=\"#000000\" points=\"8,-8 8,-83 78,-83 78,-8 8,-8\"/>\n",
       "<text text-anchor=\"middle\" x=\"43\" y=\"-15.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">n=1,...,N</text>\n",
       "</g>\n",
       "<!-- theta -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>theta</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"43\" cy=\"-129\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"39\" y=\"-125.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">θ</text>\n",
       "</g>\n",
       "<!-- Xn -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>Xn</title>\n",
       "<ellipse fill=\"#d3d3d3\" stroke=\"#000000\" cx=\"43\" cy=\"-57\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"34.5\" y=\"-54.3\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">X</text>\n",
       "<text text-anchor=\"start\" x=\"45.5\" y=\"-54.3\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\" fill=\"#000000\">n</text>\n",
       "</g>\n",
       "<!-- theta&#45;&gt;Xn -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>theta&#45;&gt;Xn</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M43,-110.8314C43,-103.131 43,-93.9743 43,-85.4166\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"46.5001,-85.4132 43,-75.4133 39.5001,-85.4133 46.5001,-85.4132\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1a18b543d0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gcp = Digraph('coin_toss_bayes_plate')\n",
    "gcp.node('theta', label='<&theta;>')\n",
    "with gcp.subgraph(name='cluster_0') as sg:\n",
    "    sg.node('Xn', label='<X<sub>n</sub>>', style='filled')\n",
    "    sg.attr(label='n=1,...,N')\n",
    "    sg.attr(labelloc='b')\n",
    "gcp.edge('theta', 'Xn')\n",
    "gcp.render('coin_toss_bayes_plate', format='png')\n",
    "gcp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To carry out Bayesian inference, we need the joint probability density of all variables.\n",
    "It is:\n",
    "$$\n",
    "p(x_1,\\dots,x_N, \\theta) = p(x_{1:N}|\\theta) p(\\theta) = \\left(\\prod_{n=1}p(x_n|\\theta)\\right)p(\\theta),\n",
    "$$\n",
    "where we first used Bayes' rule and then the fact that the coin tosses are independent.\n",
    "We must find the mathematical form of all these expressions.\n",
    "For $\\theta$, we simply have:\n",
    "$$\n",
    "p(\\theta) = 1_{[0,1]}(\\theta),\n",
    "$$\n",
    "where $1_A(x)$ is the indicator function of $A$, i.e., $1_A(x) = 1$ if $x$ is in $A$ and zero otherwise.\n",
    "For the Bernoulli pmf's we have:\n",
    "$$\n",
    "p(X_n = 1|\\theta) = \\theta,\n",
    "$$\n",
    "and\n",
    "$$\n",
    "p(X_n = 0|\\theta) = 1- \\theta.\n",
    "$$\n",
    "So, in a unified way, we can write:\n",
    "$$\n",
    "p(x_n | \\theta) = \\theta^{x_n}(1-\\theta)^{1-x_n}.\n",
    "$$\n",
    "Now, let's re-write the joint pmf:\n",
    "$$\n",
    "p(x_1,\\dots,x_N, \\theta) = \\theta^{\\sum_{n=1}^Nx_n}(1-\\theta)^{N-\\sum_{n=1}^Nx_n}1_{[0,1]}(\\theta),\n",
    "$$\n",
    "which has a nice interpretation as it depends only on the total number of heads $\\sum_{n=1}^Nx_n$.\n",
    "\n",
    "Now, we are in a position to apply Bayes rule to condition on the data.\n",
    "We have:\n",
    "$$\n",
    "p(\\theta|x_{1:N}) = \\frac{p(x_{1:N},\\theta)}{p(x_{1:N})} \\propto p(x_{1:N}, \\theta) = \\theta^{\\sum_{n=1}^Nx_n}(1-\\theta)^{N-\\sum_{n=1}^Nx_n}1_{[0,1]}(\\theta).\n",
    "$$\n",
    "It may be the first time you encounter this, but we have actually discovered a new distribution called the [Beta distribution](https://en.wikipedia.org/wiki/Beta_distribution).\n",
    "This is what the posterior turns out to be.\n",
    "This is one of the few instances where the posterior is analytically available.\n",
    "We will return to the example, once we introduce the Beta distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: The Beta distribution\n",
    "\n",
    "The Beta distribution is suitable for random variables that take values in $[0,1]$ but are not necessarily uniform.\n",
    "We write:\n",
    "$$\n",
    "X\\sim \\operatorname{Beta}(\\alpha, \\beta),\n",
    "$$\n",
    "where $\\alpha$ and $\\beta$ are positive shape parameters.\n",
    "The interpretation of the parameters is more or less this:\n",
    "+ The bigger $\\alpha$ is, the more the distribution is pulled towards zero.\n",
    "+ The bigger $\\beta$ is, the more the distribution is pulled towards one.\n",
    "\n",
    "The PDF of the Beta is:\n",
    "$$\n",
    "p(x) = \\frac{x^{\\alpha-1}(1-x)^{\\beta - 1}}{B(\\alpha,\\beta)},\n",
    "$$\n",
    "where\n",
    "$$\n",
    "B(\\alpha,\\beta) = \\frac{\\Gamma(\\alpha)\\Gamma(\\beta)}{\\Gamma(\\alpha+\\beta)},\n",
    "$$\n",
    "where $\\Gamma$ is the [Gamma function](https://en.wikipedia.org/wiki/Gamma_function).\n",
    "It's expectation is:\n",
    "$$\n",
    "\\mathbb{E}[X] = \\frac{\\alpha}{\\alpha + \\beta}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Inferring the probability of a coin toss from data (2/2)\n",
    "\n",
    "Now, that we know about the Beta distribution, we can write for the posterior of $\\theta$:\n",
    "$$\n",
    "p(\\theta|x_{1:N}) = \\operatorname{Beta}\\left(\\theta\\middle|1 + \\sum_{n=1}^Nx_n, 1 + N - \\sum_{n=1}^Nx_n\\right).\n",
    "$$\n",
    "where with $\\operatorname{Beta}(\\theta|\\alpha,\\beta)$ we mean the PDF of the $\\operatorname{Beta}(\\alpha,\\beta)$ evaluated at $\\theta$ (this is a very useful notation).\n",
    "So, we see that the $\\alpha$ parameter is just one plus the number of heads and the $\\beta$ parameter is one plus the number of tails.\n",
    "\n",
    "Let's try this out with some fake data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1a192312d0>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEtCAYAAAAPwAulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxcdb3/8dcnaZq1SbpvgAWKtLRCC8qlyCJQWhEqrUVbtkIVFUHkXhXQy5UCilbZwR8g6JUKCAJSW/YK3ELVspQihVLZu9u0aUiaNkmzfX9/nDnJTDJJZiazZt7Px2MeZ3rW70nSz3zn812OOecQEZHskZPqAoiISHIp8IuIZBkFfhGRLKPALyKSZRT4RUSyjAK/iEiWUeAXEckyCvyS0cys1MxazcyZ2V1d7POemW1McDk+bWbXmtnLZrbDzGrN7J9mdqWZFXdxTI6Z/ZeZ/cvMGsxsk5ndGG7/aPYV6YkCv2S6wwEDWoBZZhbyN21mA4CxwOoEl+PrwH8BHwLXApcB7wI/A/5hZoVhjrkZuAl4B7gEeAT4HvB4x/uIcl+RbvVLdQFEeunwwPIBYB5wDPBSh+0GvJ7gcjwK/MI5VxO07i4zex+4EvgG8Gt/g5lNwAvgjznnZget/xi4DZgL/DHafUUioZqCZLojAstrgCZgdoft/gdDQmv8zrlVHYK+70+B5cQO68/E+0C6pcP6e4A64JwY9xXpkWr8kukOBzY45z4ys+eBr5jZf7r2Saj8D4ZONf5AimRQFNeqcs61Rlm+fQLLig7rPwe0Aq8Gr3TONZjZPwPbY9lXpEeq8UvGMrMS4NO01+Yfwwu0RwbtdgTwb+fctjCn2A/YEcVrvyjLlwtcBTTTORUzCqh0zu0Nc+gWYIiZ9Y9hX5EeqcYvmWwSXuXFD/x/Ae7ES/e8Eujx8mngqS6O3wacHMX1wn14dOcW4Cjgv51z73bYVgSEC+QADUH7NEa5r0iPFPglk/lpnNUAzrkdZrYCL/BfDkwm9IMhhHOuAXguEQUzs58C3wXuds79IswudcCwLg4vCNon2n1FeqTAL5ksXMPtn4HbzWwS3eT3oS0VMzSK6+1wzrX0tJOZXQ38D/B74MIudtsKHGJm+WFSOKPxUjuNMewr0iPl+CWThcvfLwYcXq2/px49+wL/juK1b08FMrMFwALgD8AFrusnHb2G9/8vuD0CMyvAS2GtinFfkR6pxi8ZKTAgahzwTPB659wWM3sFL/C3Atudc5u7OE1cc/xmdhVwNXAfML+HHkB/Av4b+E9gRdD6b+Ll6x+IcV+RHinwS6aaBOQSvjb/Z+B6vJr/s12dIJ45fjO7GG8swcbAOc8ys+BdKpxzfw269ltm9v+A75rZY3gN0OPxRuO+SFAvoGj2FYmEAr9kqu7SOH7gty62J4Lfl34/YFGY7S8Cf+2w7j+B9cC3gFOBSuB24Kow3xai2VekW6aHrYuIZBc17oqIZBkFfhGRLKPALyKSZRT4RUSyTNr36jGzZrwPqF2pLouISIYoBVqdc2FjfES9esxsH7wnCh2B13+6GDjBObe8w37rgU+FOcUvnXM/iqrY7edsBaysrCyWw0VEsk5NTQ2Ac86FzepEWuMfi/cwiNXA88CXu9n3JeCKDuu2RHidcHaVlZWVVVdX9+IUIiLZo7y8nJqami6zJJEG/pecc8MAzGwm3Qf+T5xzL0dRRhERSaKIGnc1MlBEpO9IRK+eE81st5k1mtlbZvYd6zBpiYiIpE68e/U8gTdF7EfAYLyHQN+B9xSk/wp3gJn1lLxXq66ISBzFNfA7577bYdViM3sA+J6Z3eKc2xDP64mISPSS0Y9/EXAW3kMkOgV+51x5dwcHvhF0W+t3zlFZWUlDQwOtrWqOkOjk5ORQUFDAkCFDUFZSskEyRu7610hIRHbOsWXLFiorK2lqakrEJaSPa2pqorKyki1btqDZajPM5lVQ09VzdqQryajxz8ML+q8l4uSVlZXU1tYyfPhwBg0alIhLSBaoqqqioqKCyspKhg6N5jG8kjKbX4ffngTln4JLVkNu2k9EkDYi/kmZ2RmBt/4DJ443syHAHufc02Z2JnA68CSwGRiE17g7E7jeObcxfsVu19DQQH5+voK+9MqgQYOorq6moaEh1UWRSK0PPIWyegN8/CKMPSm15ckg0XxEPtLh31cHlhuAMcDHwBDgV3g9evYCbwHnO+fCPZEoLlpbW8nNzU3U6SWL5Obmqo0ok1SsbX//1qMK/FGIOPA757pt9QqM1p3a6xKJiEQiOPCvexxOuwnyClNXngyiaZlFJPM0N0Llu+3/bqyF95elrjwZRoFfRDJP5XvQ2uy9H/4Zb/lWx2y0dEWBvwtbqutZ9I/13P78+yz6x3q2Vtcn9fpmFtFr/fr1SS1Xshx11FF88YtfTHUxJF1tf8dbFpTD0YFxo+8tg4aa1JUpg6j/UwfbaxtYsGQtz67dRmtQl+5rHl/L9AkjuOb0CQwbUJDwcqxcuTLk31dccQXvvfceixcvDlk/cuTIhJdFJO1UvO0th0+EcadCv0Jorod1T8Dks1NbtgygwB9ke20DZ9y5ko1VdZ22tTp4+u1trN26iz9/52iGDshPaFmOOuqokH8PHDiQ/Pz8Tuu7snfvXvLzE1tGkZTxG3aHT4D8AXDwF2HtYi/do8DfI6V6gixYsjZs0A+2saqOq5a8naQSReauu+7CzFi+fDlz585l4MCBTJo0CYC5c+cybty4Tsf86Ec/oqAg9JtLS0sL119/PRMmTKCgoIChQ4cyf/58duzY0e31f/GLX9CvXz+2bt3aads555zDsGHD2kZV33///Zx00kkMHz6coqIiJk6cyMKFC3scde3f47Zt20LWP/PMM5gZL78c+giIJUuWcMwxx1BSUkJJSQnTpk3jjTfe6PYakkGCAz/AxNne8uOXvIZf6ZYCf8CW6nqeXbut5x2BZ9duS3rOPxLnnnsuo0eP5uGHH+b666+P+vizzz6bBQsWcMYZZ/D444/zy1/+kmeffZYTTzyx24FN8+bNwznHfffdF7K+traWxYsXc/bZZ5OXlwfABx98wKxZs1i0aBFPPvkk3/jGN7j++uu56KKLoi5vV+68805mzpzJgQceyEMPPcR9991HbW0txx13HO+++27PJ5D0tmcn1P7bez98orccOt5buhbY031FRZTqafPcOxUhOf3utDp4bl0F86aMSWiZojVz5kxuvPHGmI594YUX+NOf/sTvfvc7vv71r7etnzBhAkcddRQPPPAA3/jGN8IeO3r0aKZOncqiRYu44or2p24+/PDD1NXVMX/+/LZ1V199ddt75xzHHnssJSUlXHTRRdx8882UlJTEVH5fTU0Nl19+Oeeeey6LFrWPG5w6dSpjx47lZz/7WacPKMkw2/3++wbDAt9mS4Km2dhdAWWjk16sTKIaf8Cu+ugmeKupS78J4WbNmhXzsU899RT9+vXjq1/9Ks3NzW2vI444giFDhvDiiy92e/z555/PunXreOWVV9rWLVq0iMmTJ3PooYe2rXv33Xc555xz2HfffcnLyyMvL49vfetbNDc38+GHH8Zcft+KFSvYvXs355xzTsh9FBYWcuyxx/Z4H5IB/DTPoP2hf7H3Pr8U+gVSl6rx90g1/oDSwryo9i8rim7/ZOhND5+Kigqam5spLS0Nu72ysrLb42fNmkV5eTn33nsv//Ef/8FHH33E3/72N2699da2fT755BOOOeYYhgwZwk9/+lPGjh1LQUEBK1as4Pvf/z719b1Pn1VUVAAwffr0sNsLCzWyM+N1zO8DmEHxMKjZ6NX4pVsK/AFTDxnONY+vjSjdk2MwdfzwxBcqSuHmki8oKGDv3r2d1ncM5EOGDKF///6sWLGCnJzOXwTLyrp/EFpBQQFf+9rXeOihh7jllltYtGgReXl5nHXWWW37LFu2jMrKSp588kmOPPLItvXB3xK6Oz/Q6V7C3QfAb3/7Ww477LBO5wl3b5Jh2gL/xND1JQr8kVLgDxhdXsj0CSN4+u2eG3inTxjBqPLMqDmOGTOGP/7xj1RVVbXNYFpfX89zzz0Xst8pp5zCLbfcws6dOznllFNiutb555/P3XffzeLFi7nvvvuYMWMGgwcPbtvufzD179+/bV1rayv/+7//G9F9AKxZs4ZPfepTbeuXLl0ast9xxx1HcXEx7733XpdtEpLBWltg+zrvfXCNH7zAD7BbqZ6eKPAHueb0CazduqvbLp37DSri2tMndrk93Zx55pn87Gc/Y86cOfzgBz9g9+7d3HLLLZ2+HUybNo25c+cyd+5cLr30UqZMmUL//v3ZvHkzL7zwAmeccQYzZszo9lpTpkzh4IMP5rLLLmPz5s3cdtttIduPPfZYBgwYwAUXXMCCBQtoaWnhjjvuYNeuXT3ex+c//3n2339/Lr30Uurr6xkwYACPPPIIq1atCtlv4MCB3HDDDVx88cVs376dL3/5ywwaNIiKigpeeeUVysrKuOqqqyL86UnaqfrYG6gF3QR+1fh7ou+9QYYNKODR70zhlIkjyOmQNckxOGXiiKQM3oqngw8+mMWLF7N9+3ZmzZrFj3/8Y+bNm8ecOXM67fvAAw9w3XXX8cQTT/CVr3yFmTNnsnDhQgoLCznkkEMiut7555/P5s2bGT58eKcpF0aOHMnSpUsxM+bMmcPFF1/MpEmTuOGGG3o8b15eHk8++SQHHHAAF1xwAeeffz4DBw7kpptu6rTvhRdeyFNPPcXmzZuZP38+06dP5/LLL2fjxo0ce+yxEd2HpCm/R09eMZSPCd1WEki/qnG3R5buj5ozs+qysrKy6urqsNs3bPAe4xv89T8etlbX89y6CmrqmigrymPq+OEZk96R2CTqb0niaOX/g2f/G4ZNgIv+Ebrt1XvgqR/C4LFwyeupKV+aKC8vp6ampqarZ5or1dOFUeWFaddPXyTr+bX5kjCPx/Rr/Lu3J688GUqpHhHJHH7DbXG4wB/I8e/dBU3pN7I+nSjwi0jm2BNB4AfV+nugwC8imWNPIKCHC/zFQYFfDbzdUuAXkcyxJzBgL1zgzy/xevuAunT2QIFfRDKDc0GNu8PC7+M3+irV0y0FfhHJDHtroTkwPXjxkPD7qGdPRBT4RSQzBOftw6V6gtcr1dMtBX4RyQyRBP620buq8XdHgV9EMoMf+PsPgLwuRtEr1RMRBX4RyQzdjdr1qXE3Igr8IpIZuhu161ONPyIK/H3MVVddxWc/+9lO6z/55BMuvfRSRo8eTWFhIccccwxr1qzptN+YMWNCnoubThYtWsQ+++xDXV3X02ZLH9bdqF2fP4iraQ/s3Z34MmUoBf4+ZNOmTdxwww389Kc/DVlfUVHBlClTWL16NXfccQePPPIIVVVVnHbaaTQ0NKSotNE755xzKC4u5le/+lWqiyKp0N2oXV9w/3418HZJgb8PufXWWxk5cmSnJ2jNmzeP4uJinn/+eU4//XROO+00brrpJjZt2sQLL7yQotJGLzc3l29/+9vcdtttcXk+r2SY7kbt+jRfT0QU+DPAtGnTOO644/j1r3/NxIkTKSgoYMyYMSEPMGlsbOT3v/89Z599dsixy5YtY9myZdx8880hjzzcf//9Afj444+7vfaWLVsoKipi/vz5Ievvv/9+cnJyQh6mngxnnXUW1dXVPPzww0m9rqQBP5B3NWoXvN4++aWh+0snfXM+/pZm2LUl1aVoVzoacmP/Ub/++us0NjayZ88efvKTn1BeXs7tt9/OZZddxsEHH8yMGTN45ZVXqKqq4vjjjw859u6772bs2LEcffTRNDc3t633H3cY/GEQzujRo/n+97/PwoULufzyyxk/fjxLly5l/vz5XH311Vx66aUx31csRowYwbhx43jqqac477zzknptSbG2HH8Xo3Z9xUO9qZk1iKtLfTPw79oCtx6a6lK0u3QNDIztqU4fffQRVVVVjB8/nhUrVlBUVAR4z6AdMWIES5YsYcaMGaxcuRKAyZMntx3b3NzMsmXLqK2tJS8vL+z5/Zp/d6644gruuecerrzySi655BLmzJnDJZdcEvLs2gULFvDII4/wr3/9iz/+8Y/MnTs3qvv84IMPOOiggzj33HP5wx/+0LZ+zJgx3H///RxzzDFt6w4//HCWL18e1fklwzU3QkPgKXzF3dT4wevZU/WhZujshlI9ae71171HyF177bVtQR+gpKSEfffdl507dwKwdetWcnNzGThwYNs+69ato7a2luuuu47XXnst5PWVr3yFnJycsD2AOhowYABXX301ixcv5rTTTuPss8/mxhtvDNnnoIMO4tZbb+XII4+M6T7XrFnDAQccwNNPP01LSwsANTU1bNy4kUMPDf0QHzZsGNu2baO1tTWma0kGqqtsf99djh/00PUI9M0af+lor5adLkpHx3zoqlWrKCgo4PTTTw9Z75xj27ZtnHzyyQDU19eTn5+PWftT4tevXw/AlClTOgX4N998kylTplBeHvaRnJ2MGzcOgKKiIu66666Q64DX4wbguuuui/zmOpRn6tSprF69mpUrV7Z1Nx0zZgylpaUh+xYUFNDS0kJjYyMFBQUxXU8yTMh0DT2ketoCv3L8XembgT+3X8yplXSzatUqhg4d2ilV8/TTT1NdXc2pp54KwJAhQ6irq2Pv3r3k5+cDtOX0Ox67fPlyPvzwQ37yk59EVIbVq1czc+ZMPv/5z/P3v/+dBx98kHPPPTfqe7nooosAuOOOOzpte/PNN5k2bRqjR4/m8ccfbwv8HWv7AFVVVRQXFyvoZxN/8FZOPygc2P2+Cvw9Uqonza1evZqKioq2xliAvXv3cuWVV/LZz36WadOmATB+/HgAPvzww7b9DjjgAADWrl0bcuwPf/hDJk2a1FZL7866deuYPn06J554Ii+++CInnHACV155ZUz9/++4446wQR9oC/IzZsxg6dKlgPdhcNhhh3Xa96OPPmLChAlRX18yWPDgrQ7fNjspVuDviQJ/Gvvggw+orq5mxIgRzJkzh2XLlvHYY49x0kknsWnTJh588MG2lMsXvvAFAF5++eW24w877DAmTZrE1VdfzSOPPMITTzzBySefzLZt23j00UfJzc3t9vrr16/n5JNP5tBDD+Whhx4iNzeXhQsXsmnTprh249y1axcbNmzg0EMPZfLkydTW1vLBBx+EDfytra28+uqrnHDCCXG7vmSASEbt+oJn6HQucWXKYAr8aWzVqlUAPPbYY5SXlzN79mzmz5/PyJEjefXVVxk7dmzbvvvssw/HHXccS5YsCTnHY489xsSJEznvvPOYP38+Bx10EK+99hoHHnhgt9f22w9GjRrFkiVL2tJHRx55JLNnz+YXv/gFlZWV3Z4jUmvWrGG//fZry+WfdtppLFmyhLVr13ZK9SxfvpyamppO4xWkj4tk1K7Pn6itucHr1imd9M0cfx/x+uuvU1ZWxuGHH86DDz7Y4/6XXHIJZ511Ftu3b2fYMO/r7v77789f//rXqK89YsQI3n///bDbHn300U7rmpqaaGlpobW1laamJhoaGujfvz85OT3XLTrm8mfMmMHFF1+MmXX6gFq0aBHHHnssn/nMZ6K8I8lo/qjd7gZv+YoGt7+vq4KCssSUKYOpxp/GVq1axRFHHNGpB01XZs+ezeTJk1Myl803v/lNCgsLWbFiBfPmzaOwsJCXXnopZJ8LL7yQCy+8sNOxb775ZkjgP+mkk9ixYwcTJ04MufePP/6Yhx56iIULFybuRiQ9+fn6nnr0ABQE9VTz+/5LCNX405RzjjfeeINvf/vbER9jZtxzzz0888wzCSxZePfeey/33ntvt/vcddddYdf/5je/Cfl3QUEBe/bs6bTfxo0buf322zn66KNjLqdkqGhy/PmlgAEO6hX4w1HgT1NmRnV19H+0hx56aNgukH3B8ccf32lKCskSbRO0RZDqycnx0jsN1dBQk9hyZSgFfgnhD/oSSRvORVfjh6DArxp/OMrxi0h6a6iG1ibvfXePXQxWGMjzK9UTlgK/iKS33cHTNURa4w8EftX4w1LgF5H0FjxPT1EEvXqgvcavHH9YCvwikt78wF9QDv26f35EmwKlerqT8YE/JyenbRpfkd5oaWmJaMCZJFm0DbvQPmhLqZ6wMv6vvF+/fjQ1NWludukVf8Rxv37q6JZ22rpyRpjmATXu9iDjA39paSmtra3s3LkTpwmZJAbOOSorK2ltbaWsTMP70059lbcsHBT5MQXK8Xcn46s3xcXFDBgwgMrKSnbt2qUam0StubmZxsZGSktLQ55yJmmiLhD4i3qYhz9YoXr1dKdPRMlRo0ZRU1PD7t27lfKRqPXv359Bgwaptp+u/Bp/8ORrPfFz/PXV3gCwCOe7yhZ9IvDn5OQwcODAkOfNikgfURdLqicQC1wLNO6G/AHxL1cGiyjHb2b7mNmtZvY3M9ttZs7MvtDFvmeZ2Ztm1mBmm81soZnpGXkiEpu2Gn8Ugb8waIZONfB2Emnj7ljgTGA38HxXO5nZOcADwN+BU4CfAxcD9/aqlCKSveo+8ZaxNO6CGnjDiDTV85JzbhiAmc0EvtxxBzPLBa4HljrnLgqs/j8zawLuNrObnXOvxKPQIpIlmhuhsdZ7H02NP/jhK2rg7SSiGr9zLpIW06OAEcCiDusfAJqA2dEVTUSynp/mgehq/Ln9oH9J4BwK/B3Fsx//xMDy7eCVzrk64MOg7SIikakLCvzR1PhBE7V1I569evy+VlVhtlUFbQ9hZj39VtTHTiRbhdT4o+y1V1gOuzYrxx9GIkbudjV8VsNqRSQ6fo0/vwxy86I7VhO1dSmeNf6dgeXgoPe+QcDH4Q5yzpWHW+8LfCNQrV8kG9XHMGrXp4nauhTPGv/awDIkl29mRcCBdMj9i4j0KJbBWz5N1NaleAb+l4FtwLkd1p8J5AGPxfFaIpINYhm85dNEbV2KONVjZmcE3n4usDzezIYAe5xzTzvnms3sR8C9ZvZr4FFgPPBL4FHn3MvxLLiIZIFYBm/5NFFbl6LJ8T/S4d9XB5YbgDEAzrlFZtYCXAF8E6gE7gIW9KqUIpKd6gLNhTHV+IMmapMQEQd+51xE09s55+4H7o+5RCIivljm4vepH3+XMv5BLCLSh9X1IsevB653SYFfRNJXPBp3mxugqSF+ZeoDFPhFJD21tkJ9HBp3QemeDhT4RSQ97a0Bf37I3jTughp4O1DgF5H0VBfjzJy+AtX4u6LALyLpyU/zQGw1/rwC6Bd4+J8aeEMo8ItIevL78OfmQ15RbOfQRG1hKfCLSHoK7sppEQ0j6kwTtYWlwC8i6amtK2fYR3lERhO1haXALyLpqW1mzhimZPZporawFPhFJD31ZvCWTxO1haXALyLpqTdz8fs0UVtYCvwikp7iUePXRG1hKfCLSHrqzVz8Pk3UFpYCv4ikp97Mxe9TP/6wFPhFJD31Zi5+nxp3w1LgF5H001jnTacMvevH7zfuNu6Glqbel6uPUOAXkfRTHzRBWzxSPQANu2I/Tx+jwC8i6SdkZs5eDODSnPxhKfCLSPrxa/yWE1prj1bwsWrgbaPALyLpx6/xF5RDTi/CVP9iyOnnvQ+e5jnLKfCLSPqJx+At8Gb19FNFSvW0UeAXkfQTj+kafG19+VXj9ynwi0j6aRu81YuunD5NzdyJAr+IpJ89ld6yOB6BX6mejhT4RST91AUCf9GQ3p9L0zZ0osAvIulnTyDVUxyHwF+oHH9HCvwikn4SUeNXqqeNAr+IpBfngnL8Q3t/Pj/Hr1RPGwV+EUkve3dBa2BCtbg07qrG35ECv4ikF7+2D3Fu3FWO36fALyLpJTjwx6VxN5DqaaqD5sben68PUOAXkfTiN+zmFUNeYe/Ppxk6O1HgF5H0Es/BW9Bhhk6le0CBX0TSTTy7ckJojV89ewAFfhFJN22Dt+LQlRO8dFG/Au+9Uj2AAr+IpBu/xh+Phl2fevaEUOAXkfTi5/jjMTOnTzN0hlDgF5H0smeHt4xnjV8zdIZQ4BeR9NI2F79SPYmiwC8i6SNknp541viV6gmmwC8i6aNxN7Ts9d7Hs8avVE8IBX4RSR/xnq7Bp4exhFDgF5H04ef3IUGpHuX4QYFfRNKJX+PvVwj9i+N3Xj2MJYQCv4ikj0R05YTQh7E4F99zZyAFfhFJH3UJGLwF7amelr3QVB/fc2cgBX4RSR+J6MoJoTN0Kt2jwC8iaSQRg7egPdUD6tmDAr+IpJNE1fgLNSd/MAV+EUkfiZiZEyA3z3uiFyjVgwK/iKSTPQlK9UBoz54sp8AvIukjUd05QYO4gijwi0h6aNwDzYGulomo8WsQV5u4Bn4z+4KZuS5e4+J5LRHpY0Lm6YlzP37QDJ1B+iXovFcAL3VYtz5B1xKRvqAuOPDH6Xm7wQpV4/clKvC/55x7OUHnFpG+yG/Yzc2H/iXxP78extJGOX4RSQ/BXTnN4n9+pXraJCrw/8bMms2sxsyeMLMjEnQdEekrEvGQ9WB6GEubeKd6aoBbgOVAFTAe+BHwdzM73jn3SscDzKyn30JZnMsoIukoUYO3fEr1tIlr4HfOvQG8EbRqhZktBd4GrgOmxvN6ItKH7Pb78CegYRdCUz3OJSadlCES1bjbxjm3zcyWAV/uYnt5uPW+wDcC1fpF+rpdW7xl6ajEnN9P9bgW79m++QMSc50MkKzG3RxATz8Qka7t2uotS0cn5vwFmqjNl/DAb2YjgJMBde8UkfCcCwr8Ca7xQ9b37IlrqsfMHgA+AlYDnwDj8AZzFQI/jue1RKQPaaiBpj3e+0QF/oIywAAH9VWJuUaGiHeO/y1gLnAJUAzsxOvh8zPn3NtxvpaI9BV+bR8Sl+rJyfV6DO3Z0d6QnKXi3atnIbAwnucUkSzgB/6cvMRM0OYrHhYI/BWJu0YG0MhdEUm9th49IyEngWGpZJi3VOAXEUmxRPfo8ZUM95Z7sjvVo8AvIqmX6D78vpLA4DDV+EVEUizRXTl9fo0/yxt3FfhFJPWSnepRjV9EJMWSVeP35wGqq4TWlsReK40p8ItIau2thb013vtk1fhdK9TtTOy10pgCv4ik1q5/t+5V/5oAABAsSURBVL9PVo4fsjrdo8AvIqnl9+ix3NDAnAiFAyEnMG5VgV9EJEX8/P6AEd60ComUk9Oe58/inj0K/CKSWslq2PVp9K4Cv4ikWLIGb/k0eleBX0RSLFl9+H3FqvEr8ItIainVk3QK/CKSWqlK9ahxV0QkBZrq25+GlaxUjyZqU+AXkRQKfvLWgJHJuaZf46+vgpam5FwzzSjwi0jqpDLwQ9b27FHgF5HU8QN/8TDo1z851/QHcEHWpnsU+EUkdZLdsAtQUAa5+d773duTd900osAvIqmT7D78AGZBPXsU+EVEkqt6o7csS2Lgh6zv2aPALyKps+Nf3nLIp5N73SyftkGBX0RSo3FPe41/6LjkXjvLR+8q8ItIalS+BzjvfbIDf9t8Pcrxi4gkz/ZAmqdocHvOPVlKFPhFRJJvxzpvOXR88q+tXj0iIing1/iHJTnNA+01/r010NSQ/OunmAK/iKSG36Mn2fl9aA/8AHuyr9avwC8iyde4B6o3eO9TEfiLgwJ/FqZ7FPhFJPl2vNv+flgKcvz5JdC/xHsfPFFcllDgF5Hk89M8RUOgeEhqyjDoAG+58/3UXD+FFPhFJPn8wJ+K2r7PTzEFf/vIEgr8IpJ821PYsOsberC39D+EsogCv4gkX1sf/oNTV4a2Gv970NqaunKkgAK/iCTX3t3tc/SkQ6qnub69h1GWUOAXkeSqfK/9fSpG7foGjoHcwFO/sizPr8AvIsnl59SLh0Lx4NSVI7cfDD4otExZQoFfRJJru5/fT2HDrq+tgVc1fhGRxNn+jrdMh8DvtzGoxi8ikiDNe2HDSu/96CNSWxYIrfE7l9qyJJECv4gkz4Z/QNMe7/3YqaktC7R/62jaAzWbUluWJFLgF5Hk+eA5bznq8OQ/fCWcQQdATj/vfRbl+RX4RSR53l/mLQ+altpy+HLzYPBY730W5fkV+EUkOao+bu/Dny6BH7Jy6gYFfhFJDj/NUzQERk1ObVmCZeFkbQr8IpIcfppn7FTISaPQk4U9e9Lopy8ifVZTPXz8kvf+oJNTW5aO/Br/3l1Z81AWBX4RSbz1f4PmBrAcOPDEVJcm1OCxXrmgfdbQPk6BX0QSz0/z7PM5KBqU2rJ01C+/vdb/4f+ltixJosAvIonVsAvefsx7n25pHt/E2d7yzQehuTG1ZUkCBX4RSay/3Qx1lZBXBJPOSXVpwpt0lpfuqdsJ7z2d6tIknAK/iCRO9SZ4+Q7v/ecvhdKRqS1PV0pHtU8h8cb9qS1LEijwi0jiPH+t16hbMgKOviTVpene5HO95QfP9fnePf1SXQAR6aO2vA5vPey9P/F/oH9xbKeprue5dyrYVd9EaWEeJx8ynFHlhXEsaMCnv+gNLqurhH8+AMddFv9rpIm4B34zKwF+DnwVKAfWAtc655bG+1oikqZ274DHL/XeD5/o5dCjtL22gQVL1vLs2m20Bo2ruubxtUyfMIJrTp/AsAEFcSow0K8/HDYXVv7aS/cc84P0GmgWR4m4q8XA2cD/AKcC7wCLzexLCbiWiKSbinfgnhNh21uAwfSfQ05uVKfYXtvAGXeu5Om3Q4M+QKuDp9/exhl3rmRH7d74lRva0z2frId1S+J77jQS18AfCO5TgQucc79zzr0AnAesBG6M57VEJM00N8Kah+F306Bmo9eLZ879cMDxUZ9qwZK1bKyq63afjVV1XLXk7VhLG96wcbDf0d77R78OL9/ZJ6dxiHeqZxZQA7R9VDrnnJktAu42s0Occ+/E+ZoikgpNDVC9AXZ+6A3QWrsYGqq9bQNGwVkPwcjDoj7tlup6nl27LaJ9n127ja3V9fHN+c/+LTx0Jvz7TXjmR943l6O+A0MO9tJBfUC8A/9E4B3nXGuH9WuCtwdvMLPqHs5ZFlNJNvzD6z8sIu26rL26Dvu49mVrC7Q2e6+methbC427oa4q9Djw+sIfNB1OuznmrpvPvVPRKb3TlVYHz62rYN6UMTFdK6yy0TD/GVj6XXj7z15D7z8f8B7YMngs5A+A3HzvQ8CCkyYW9NY6nTZmJ/x33GczjXfgHwy8F2Z9VdD25Nhd0T5MXEQSJzcfRkyEiWd4I2AHDO/V6XbVN0W1f01ddPtHpH8RzP6d943l77d6A7tam1MzZ/+R3477KRPRnbO7z+pO25xz5d2dLPCNIPpa/+Cx6d9vWCQluqiNhtRSLfBv82q6Of28Btq8QuhfAvklXtfHwQd6aZ049n4pLcyLav+youj2j5iZN+js6O9B7b+hYi1Uvg9NddDS6D043g9pId+k4twmUL5ffM9H/AP/TsLX6v1ZmarCbEuMEZ/xXiKSUaYeMpxrHl8bUbonx2Dq+N59w+iRmTeyt3RU+s41FKV4d+dcC4w3s47n9SNwnJvgRaSvGV1eyPQJIyLad/qEEYkZzNXHxTvwL8YbtDWjw/p5wLvq0SMikbjm9AnsN6io2332G1TEtadPTFKJ+pZ4B/6ngP8DfmdmXzezE8zsXuAYoO+OfxaRuBo2oIBHvzOFUyaOIKdDk0SOwSkTR/Dn7xzN0AH5qSlghjMX58EJZlaKN2XDGXi1/3fwpmz4S4znqy4rKyurru6p16eI9EVbq+t5bl0FNXVNlBXlMXV8gubq6UPKy8upqamp6arzTNwDf7wp8IuIRKenwN83ZyASEZEuZUKNvxWwsrLYBvCKiGSbmpoa8GbMCVu5z4TA34z3zWRXDIf7nxY18StR2tM9933Zdr+ge45WKdDqnAs7VivtA39v+PMA9TQ6uC/RPfd92Xa/oHuO97mV4xcRyTIK/CIiWUaBX0Qkyyjwi4hkGQV+EZEso8AvIpJlFPhFRLJMn+7HLyIinanGLyKSZRT4RUSyjAK/iEiWycjAb2YlZnabmf3bzOrNbJWZfTnCYw80s7+YWY2Z1ZrZU2Z2SKLL3Fux3rOZXWBmS81sQ+C49wPnGZqMcvdGb37PQecwM3vBzJyZ3ZKossZDL/+uzcy+ZWavm1mdmVWb2ctmdnSiy90bvbzn2Wb2DzP7JPBaaWZfS3SZe8vM9jGzW83sb2a2O/C3+YUojj/CzJ43sz2B+37IzEZHU4aMDPx4z/Y9G/gf4FS8p3wtNrMvdXeQmQ0DVgBjgPOAM4FBwItmtk8iCxwHMd0zcA3ezKY/Br4I3AR8DXjNzNJ9wqtY7znYN4FxCShbIvTmfn8L/Ar4M/ClwHmeAooTU9S4ifX/8nnAo8BW4KzAawvwJzP7ekJL3Htj8WLPbuD5aA40s/HAcsDwnnL4TWAysNzMSiI+kXMuo154f9QOmBW0zoC/Aet6OPZXQD0wKmjdYLzAeGeq7y1B9zwszLrjA+e7JNX3loh7Dtp/NFANzA6c65ZU31eCfsezgRZgSqrvI4n3vBxYD+QErcsJrFue6nvroezBZZ4Z+Bl8IcJjH8b7sCsOWjcu8Pu/ItIyZGKNfxbe/NRL/BXOu/tFwLge0jazgL8657YGHbsTeBz4SmKKGxcx37NzbnuY1a8Flun8Lac3v2ffncBLzrk/J6aIcdWb+70E7z5XJraIcdebe24CdjvnWoOObcWrRe9NTHHjI7jM0TCzPOA04FHn3J6g8/0LeBmvAhCRTAz8E4F3wvzw1gRt78TMCoEDgbfDbF4DDAukgtJRTPfcjRMDy3A/i3TRq3s2szOBE4CLE1C2RIj17zoPOAp4y8x+bmYVZtZsZmsD6ZB01pvf8a+B8WZ2pZkNMbOhZnYlcDBwcwLKmg4OAArpOoZFHAcyMfAPBqrCrK8K2h7OQLyvkbEcm2qx3nMnZjYIuA14H+9rY7qK+Z7NbAhwK3Clc25TAsqWCLHe72AgH6/N6nTgu8ApwFvAvWb2zTiXM55i/h0755YAXwZ+COwAtuO1Y33VOfdMnMuZLvyfR1c/s8JABbdHYR/LlQG6G27c01Dk3hybSr0ut5kVAX/Ba9A+zjmX1l+Jif2ebwM+xqsVZpJY7tevvBUAX3LObQAws+fwaohXAffErYTxF9Pv2MxOBv4IPIjXoJ2L10j8oJmd4Zx7Mq6lTC+9jgWZGPh3Er4mMCiwDPdpCPAJ3g8llmNTLdZ7bhOoCSzF6wEw3Tm3podDUi2mew4EhDl46axSMwvenB/oybTbOdccx7LGQ2//rv/lB33wcuVm9gzwEzMb1kVbT6rF+js2vHaAF5xzFwZteibQO+92oC8G/p2BZVc/s3rnXEMkJ8rEVM9avNxex7J/JrAMm7d2ztUDHxE+D/YZYEea/ueAGO/ZZ2YFeA1oU4DTnHP/iH8R4y7We56A93e9HC8o+i+ACwPvp8a1pPHRm7/rD7o4p/+pF1NjYhLE+jseDowEVoXZtgrYP/A339d8hNcrsasYFnGbXSYG/sVAOTCjw/p5wLvOuXd6OPZkMxvhrwjkvGcAj8W7oHEU8z2bWT5eeudY4HTn3IsJK2V8xXrPj+I16nZ8gZcSOAF4Ne6l7b3e/F0/hhdAx/grArXiU4CPnHOV8S1q3MR6z58ADcCRYbYdBeyMtOabSZxzTXjfZGYH0rYAmNmn8Sp1kcewVPdpjaEPrAEvAJXA1/H+I9+LV6uZEbTfcgK9w4LWDQe2AavxGsJOBVbifYXaL9X3lqB7fhwvFXAN3n+K4NeBqb63RNxzF+dL9378vfkdDwY2Af/CGxh0Ct4HoAPmpPreEnTPNwfu77d4AxNPBf4UWHdlqu8tgns/I/D6ZaDMCwL/PiVon/XA+g7HHYLXZfW5wH3PDvzePwQGRHz9VP8AYvyhleI13G3D++RfDczssE/YgAAchJf22BX4AT4NTEj1PSXqngN/VF297k31fSXq9xzmXGkd+Ht7v3ij0R+hvTb8Wsdj0/HVi7/rXODbwOt4g/Sq8Pqyn0Nguvl0fnXzf3J90D6dAn9g/efwPjD3BO79YWDfaK6v+fhFRLJMJub4RUSkFxT4RUSyjAK/iEiWUeAXEckyCvwiIllGgV9EJMso8IvEyMymm9nywOPzdpjZr/voVAHSxyjwi8TAzH4APAP8G/gvvBHSF+NNBy2S1jSASyRKZjYVWAZc7py7IWj9M3jTDgx1zu1KVflEeqIav0gUAjNJ3gq8AdzYYfNyoD/RPxFNJKkycT5+kVSajjdR1vmu89flxsCyLLlFEomOAr9IdOYALcCKwCMegw0PLGuTWySR6CjHLxIFM9sA7NfDbqOdc1uTUR6RWCjwi0QoUMPfgfcAkTvC7PIwsNc5NzKpBROJklI9IpE7ILB8zTn3XPAGM9sfGIj3AHCRtKZePSKRKwksw+Xwzwgs/5SksojETIFfJHJ+3/zS4JVm1h/4DvAu3jNRRdKaAr9I5N4B6vC6dAa7Du/Rh99zzrUku1Ai0VKOXyRCzrk6M/st8D0zux94Ee/B5rOAy5xzy1JaQJEIqVePSBQCaZ1fAWcDRXgP+/65c+6ZlBZMJAoK/CIiWUY5fhGRLKPALyKSZRT4RUSyjAK/iEiWUeAXEckyCvwiIllGgV9EJMso8IuIZBkFfhGRLKPALyKSZf4/7lKrauy77mQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Take a fake coin which is a little bit biased\n",
    "theta_true = 0.8\n",
    "# This is the random variable corresponding to a coin toss\n",
    "X = st.bernoulli(theta_true)\n",
    "# Sample from it a number of times to generate our data = (x1, ..., xN)\n",
    "N = 200\n",
    "data = X.rvs(size=N)\n",
    "# Now we are ready to calculate the posterior which the Beta we have above\n",
    "alpha = 1.0 + data.sum()\n",
    "beta = 1.0 + N - data.sum()\n",
    "Theta_post = st.beta(alpha, beta)\n",
    "# Now we can plot the posterior PDF for theta\n",
    "fig, ax = plt.subplots()\n",
    "thetas = np.linspace(0, 1, 100)\n",
    "ax.plot([theta_true], [0.0], 'o', markeredgewidth=2, markersize=10, label='True value')\n",
    "ax.plot(thetas, Theta_post.pdf(thetas), label=r'$p(\\theta|x_{1:N})$')\n",
    "ax.set_xlabel(r'$\\theta$')\n",
    "ax.set_title('$N={0:d}$'.format(N))\n",
    "plt.legend(loc='best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions\n",
    "\n",
    "+ Experiment with different values of $\\theta_{\\text{true}}$ and different values of $N$.\n",
    "+ Is the true value always covered by the posterior PDF?"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  },
  "widgets": {
   "state": {
    "968a18c908c6495081a9936faebe0473": {
     "views": [
      {
       "cell_index": 11
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
